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We develop a variational scheme called “Gutzwiller renormalization group” (GRG), which enables 
us to calculate the ground state of Anderson impurity models (AIM) with arbitrary numerical 
precision. Our method can exploit the low-entanglement property of the ground state in combination 
with the framework of the Gutzwiller wavefunction, and suggests that the ground state of the AIM 
has a very simple structure, which can be represented very accurately in terms of a surprisingly small 
number of variational parameters. We perform benchmark calculations of the single-band AIM that 
validate our theory and indicate that the GRG might enable us to study complex systems beyond 
the reach of the other methods presently available and pave the way to interesting generalizations, 
e.g., to nonequilibrium transport in nanostructures. 
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Introduction .— Impurity models are ubiquitous in con¬ 
densed matter theory. Originally, the AIM was deviced to 
describe magnetic impurities embedded in metallic hosts 
and heavy-fermion compound mm- However, the AIM 
can be applied to describe many other physical systems, 
such as quantum dots mg and dissipative two-level sys¬ 
tems [g. The importance of impurity models in con¬ 
densed matter has grown further with the emergence of 
dynamical mean-field theory (DMFT) [7] and its success 
in describing strongly correlated materials mil]. In fact, 
within DMFT, solving the many-body lattice problem 
amounts to solve recursively an auxiliary AIM. 

Among the many methodologies developed to solve im¬ 
purity models, one of the most successful is the numerical 
renormalization group (NRG) |21 [H], which was origi¬ 
nally deviced by Wilson in order to study the Kondo 
problem, and was based on the idea of “scaling”, previ¬ 
ously introduced by Anderson. The starting point of the 
NRG procedure consists in representing the AIM as a 
one-dimensional semi-infinite tight-binding linear chain 
with the impurity situated at one of the boundaries. 
Within this representation, and by using a logarithmic 
discretization of the density of states corresponding to ex¬ 
ponentially decaying hoppings (Wilson chains), the NRG 
algorithm consists in a truncation scheme that enables 
to take into account the coupling between the different 
length scales progressively, starting from near the im¬ 
purity (high energies) to longer distances (low energies), 
while retaining a relatively small number of states at each 
stage. 

More recently, the NRG algorithm has been reinter¬ 
preted as a variational method within the framework of 
the “matrix product states” mm]. From this perspec¬ 
tive, the accuracy of NRG can be traced back into the 
fact that, due to the locality of the interactions within the 
above mentioned one-dimensional representation of the 
AIM, the ground state has very low entanglement, and 


can be consequently accurately represented by a matrix 
product state with a small bond dimension. Thanks to 
this interpretation it was also possible to uncover a con¬ 
nection between NRG and the density matrix renormal¬ 
ization group (DMRG) method [HI m j which can also 
be formulated as a variational approximation within the 
framework of the matrix product states m- Eventually, 
these ideas resulted in an entirely new class of methods 
named “variational renormalization group” (VRG) |I8j . 
with broader applications with respect to both NRG 
and DMRG, encompassing, e.g., systems with dimension 
higher than I. 

In this work we introduce a variational technique 
to solve the AIM called “Gutzwiller renormalization 
group” (GRG), that combines the Gutzwiller variational 
method pRtiST] with some of the above mentioned ideas 
underlying NRG and the other VRG methods. Similarly 
to the ordinary Gutzwiller theory, the GRG consists in 
optimizing variationally a wavefunction represented as 
a Slater determinant extended over the entire system — 
which is infinite — multiplied by a local operator (named 
“Gutzwiller projector”), which modifies the correspond¬ 
ing electron configurations. However, in GRG the action 
of the Gutzwiller projector is not limited only to the cor¬ 
related impurity. In fact, in analogy with the NRG nu¬ 
merical procedure, in GRG the longer length scales are 
taken into account with arbitrary accuracy by extending 
progressively the action of the Gutzwiller projector also 
to a portion of the one-dimensional bath lying nearby 
the impurity. Thus, the size of the “projected” region is 
increased until the desired level of precision is reached. 

We perform numerical calculations that demonstrate 
the quality of the GRG variational ansatz, and present 
arguments concerning the scaling of the computational 
complexity of our algorithm indicating that the GRG 
might enable us to solve AIM beyond the reach of any 
other presently available technique. Another appealing 
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Figure 1: Upper part. Representation of a single-band An¬ 
derson impurity model with first nearest neighbor hopping. 
Each site (bath and impurity) has the elementary spin degree 
of freedom a — ±1/2. Lower part. Equivalent representa¬ 
tion of the same model where the impurity and a portion of 
the bath are formally considered as a larger impurity with 
effective elementary degree of freedom E. 


property of our approach is that it enables us to describe 
the ground state of the AIM in the thermodynamical 
limit and at all length scales, while this is technically 
very difficult with any other existing technique, includ¬ 
ing NRG, DMRG, and the Bethe ansatz P^24j . 

The Anderson Impurity Model .— For sake of simplic¬ 
ity, here our approach will be formulated for the single¬ 
band AIM represented in Fig. where the correlated 
impurity is connected to a bath constituted by a linear 
chain with first nearest-neighbor hopping: 
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The generalization of our theory to generic impurity mod¬ 
els is straightforward. 

For later convenience, we point out that the hybridiza¬ 
tion function of Himp with respect to the correlated im¬ 
purity is given by: 
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D = 2t, r = and sgn(w) is the sign function. 

The Method .— In order to illustrate our method, it 
is convenient to observe that the subsystem of I?imp 
consisting of the impurity and the first nearest N bath 
sites can be equivalently regarded as a larger impurity 


Oj.}], where the label S = 1,.., 2(A±1) runs over 
the sites 0,.., N and the corresponding spins cr = ±|, see 
Fig. 0 Within this definition, iLimp can be schematically 
represented as follows: 

Himp — ®i;}] ^ 

bb' 

+ ([Twit'S > (4) 

bH 

where the hybridization function Ajy of Sn 

^n{z) ^ Vj, vl (5) 

z — tSN 

is a well defined M x M matrix, where M = 2{N ± 1). 

Our approach consists in optimizing variationally a 
Gutzwiller wavefunction represented as follows: 


|^^) = iP^|4-o), (6) 

where |4'o) is a generic Slater determinant, and "Pat, 
which is called “Gutzwiller projector”, is the most gen¬ 
eral operator acting within the subsystem Sn- As in 
the ordinary Gutzwiller approximation scheme, in order 
to simplify the task of calculating the expectation value 
of the Hamiltonian [Eq. , the variational freedom is 
further restricted by assuming the following conditions: 

(4'o|lPi,lPjvl^o) = (^'o|^o) = l (7) 

ala^, Ivko) = ('kol a^a^, |vI/o) V E, E', (8) 

which are called “Gutzwiller constraints”. We point out 
that, while for bulk systems the Gutzwiller variational 
ansatz is supplemented by the so called “Gutzwiller 
approximation” |21j . for impurity models the method 
is purely variational, and no further approximation is 
needed 

For N = 0 the procedure described above reduces to 
applying the Gutzwiller projector only onto the corre¬ 
lated impurity — which is the standard Gutzwiller ap¬ 
proach for the AIM used, e.g., in Refs. [26l[27]. In this 
work, instead, the region of action of the Gutzwiller pro¬ 
jector is extended systematically by increasing N until 
the desired level of accuracy is obtained. 

The most evident difficulty to be overcome is that 
the number of complex independent parameters defin¬ 
ing Vn scales as 2^'^. In fact, this makes the task of 
minimizing the total energy with respect to the wave- 
function [Eq. ([^] nontrivial already for small N. How¬ 
ever, fortunately, this technical problem can be efficiently 
solved even for relatively large N thanks to the method 
of Refs. [25l |28l [29], which is summarized in the supple¬ 
mental material for completeness. A remarkable aspect 
of this numerical scheme is that it enables us to map the 
above nonlinear constrained minimization into the much 
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simpler task of calculating iteratively the ground state 
1$) of a finite AIM represented as follows: 

= SN[{ala^}] 

M M 

+ E +H.C.) + ^ (9) 
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at half-filling, i.e., within the subspace such that: 


r M 


-5:^1 


M 


T.fsfs 


1$) = M|$). 


( 10 ) 


The complex coefficients Pgs and AJ^,, have to be deter¬ 
mined numerically following the procedure summarized 
in the supplemental material. The number of AIM rep¬ 
resented as in Eq. (|^ to be solved in order to converge 
scales as o(M^), and they can be solved independently 
(in parallel) if necessary. 

As shown in the supplemental material, the ground 
state 1$) of converged parameters V and 

A'^ encodes the expectation value with respect to the cor¬ 
responding Gutzwiller wavefunction Eq. ([^ of any ob¬ 
servable A in Sn, i.e.: 

(^'jvl i[{a^, a^}] I'Sn) = (^>1 a^}] |$). (11) 

For this reason, as in Ref. [53] , here we call “embed¬ 
ding Hamiltonian”. Concerning the fact that in Eq. 
the subsystem Sn is coupled with a bath of equal dimen¬ 
sion, it is interesting to observe that, at least in principle, 
a Hamiltonian whose bath has the same dimension of the 
impurity is sufficient to represent exactly the impurity 
ground-state properties of any system (no matter how 
big it is). This fact can be readily demonstrated making 
use of the Schmidt decomposition [sniisi]. 

Let us now discuss how the computational effort to cal¬ 
culate the ground state of scales with M. In order 

to answer this question it is important to note that iJ®™^ 
is quadratic for all degrees of freedom made exception 
for those corresponding to {cl ^,Cq within the original 
representation of the AIM, see Eq. Q. Thanks to this 
observation, Eq. ([^ can be always transformed into a fi¬ 
nite one-dimensional chain by tridiagonalization [5] . This 
enables us to calculate |<i)) using efficient techniques such 
as DMRG [T31 [Ml , whose computational cost grows only 
polynomially with M rather than exponentially. Note 
that in order to perform the specific calculations that we 
are going to discuss in this work it has not been actu¬ 
ally necessary to resort to this stratagem. However, this 
might be needed in order to study impurity models more 
complicated than Eq. Q. 

We point out that the possibility to incorporate DMRG 
within our algorithm constitutes a further connection be¬ 
tween GRG and the VRG methods mentioned in the 
introduction (see also Refs. [311 I33|), as it reflects the 
fact that our approach can enforce the low-entanglement 



Figure 2: (Color online) Convergence of the GRG expectation 
value of the double occupancy d with respect to the number of 
conduction-bath sites N included within the projected region. 
The result is shown for different U/T and DjU in comparison 
with CTQMC (black crosses). 


property of the ground state within the framework of 
the Gutzwiller wavefunction, and exploit simultaneously 
both of these ideas [53] . 

From the considerations above we deduce that the 
computational complexity of the GRG algorithm scales 
polynomially with N, i.e., with the size of the region of 
action of the Gutzwiller projector. Our remaining task is 
to understand how big N has to be in order to describe 
accurately the ground state of the AIM. 

Benchmark calculations .— In order to assess the qual¬ 
ity of the GRG variational ansatz, here we perform 
benchmark calculations of the single-band AIM at half¬ 
filling, see Eq. Q. 

In Fig. is shown the behavior of the impurity dou¬ 
ble occupancy d as a function oi D/U for different val¬ 
ues of U/T. The GRG results obtained with different 
values of N are shown in comparison with continuous 
time quantum Monte Carlo (CTQMC) [34| [35] as imple¬ 
mented in TRIQS [36] . which were obtained at tempera¬ 
ture T/F = 0.02. From these calculations it emerges that 
the A^ = 4 GRG wavefunction is already sufficient to re- 
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Figure 3: (Color online) Upper panel: Convergence as a func¬ 
tion of N of the GRG double occupancy d for D/U — 4 and 
U/V = 3 in comparison with exact diagonalization for the 
truncated system composed by 2{N -|- 1) sites. Lower panel: 
Corresponding convergence of the total energy £. The results 
calculated with CTQMC are indicated by the dotted horizon¬ 
tal lines. 


produce very accurately the CTQMC double occupancy 
for all of the interaction parameters considered. How¬ 
ever, we observe that the convergence of d with respect 
to N becomes increasingly slower for larger D/U, while 
it is not very sensitive to U/T. It is interesting to com¬ 
pare this trend with the behavior of the dimension of 
the Kondo cloud. From the Bethe ansatz solution of the 
single band AIM we know that, in the large-13 limit, 
is given by [2]: 


- _ _ D 2U + 

= Vf/Tk ~ Y 1^®*^ "" ^ ^ 


( 12 ) 


(in units of the lattice spacing). As we see from Eq. (12), 


diverges exponentially with U/T. Furthermore, we 
note that is much bigger than the value of N neces¬ 
sary to achieve convergence within the GRG. In fact, for 
instance, substituting the values D/U = IQ and U/T = b 
in Eq. (12) we find that ^ 275 » 5. Thus, we deduce 


that these two length scales are not directly related, as 
one might naively expect. 

In Eig. [^is shown the behavior of the impurity double 
occupancy and of the total energy calculated with GRG 
as a function of the convergence parameter N. Note that 
the zero of the total energy is conventionally assumed 
to be the ground state energy of Eq. Q with t' = 0. 
These data are shown in comparison with GTQMG and 
with the exact diagonalization results obtained by taking 
into account only the first 2{N -\- 1) sites of Eq. 0 (and 
discarding the others). 


We observe that both the GRG energy and double oc¬ 
cupancy converge very rapidly to their respective exact 
values. From this observation we argue that the above 
mentioned GRG convergence concerns the entire ground 
state of the AIM, and not only the expectation value of 
the impurity degrees of freedom. 

From the comparison between GRG and exact diag¬ 
onalization it emerges that the dimension of the aux¬ 
iliary AIM [Eq. (|^] to be solved in order to calculate 
the GRG solution of any order N is much smaller than 
the dimension of the truncated AIM that would enable 
us to solve the problem with comparable precision di¬ 
rectly with DMRG [T3]. Thus, as discussed before, using 
DMRG to solve iteratively Eq. ([^ within the GRG algo¬ 
rithm seems to be a more convenient option. 

Conclusions .— In summary, we have developed a vari¬ 
ational method called GRG, that takes advantage si¬ 
multaneously of two among the most powerful ideas in 
condensed matter theory: (1) the Gutzwiller wavefunc- 
tion, that enables us to incorporate the most general 
single-particle wavefunction (Slater determinant) within 
the variational space, thus reducing the many-body prob¬ 
lem to correcting variationally the corresponding electron 
configurations with the “Gutzwiller projector”; (2) the 
VRG methods, that enable us to exploit the fact that 
the ground state of the AIM has low-entanglement (once 
the problem is formulated in one dimension). Using the 
GRG, we have shown that the ground state of the AIM 
has a very simple structure, which can be represented 
very accurately in terms of a surprisingly small number 
of variational parameters. These insights resulted in an 
efficient algorithm that might enable us to study complex 
systems beyond the reach of any other method presently 
available. Another remarkable property of our approach 
is that it enables us to describe the ground state of the 
AIM directly in the thermodynamical limit and at all 
length scales, while this is technically very difficult with 
the other VRG methods, which require to truncate part 
of the chain. Our work paves the way to several gen¬ 
eralizations. In particular, it will be very interesting to 
generalize it to finite temperatures m and to nonequi¬ 
librium transport in nanostructures, see Refs. HSIIIII- 
In fact, the physics of the AIM out of equilibrium is not 
still well understood, and different methods seem to pro¬ 
duce different results even for simple systems such as the 
single band AIM [35]. 
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